{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# M3.D4: Thermodynamic potentials for spontaneity and other thermodynamic properties\n", "\n", "## Learning Objectives\n", "\n", "* Calculate thermodynamic potentials using PySCF of several diatomic ideal gases\n", "\n", "## Ensembles and variables and Thermodynamic potentials for the 2nd law\n", "\n", "Thermodynamic functions can be calculated as a function of its partition function. We will just summarize the most relevant thermodynamic potentials down here.\n", "\n", "| Extensive property | Corresponding intensive property |\n", "|------------------------------|------------------------------------------------------------|\n", "| (depends on the size of the system) | (it can be measured at any point of the system) |\n", "| Volume (V) | Pressure (p) |\n", "| Energy (E) | Temperature (T) |\n", "| Moles (N) | Chemical potential ($\\mu$) |\n", "\n", " \n", "\n", "| Ensembles | Potential measuring spontaneity |\n", "|-----------------------------|-------------------------------------------------------------|\n", "| Microcanonical (NVE) | Entropy $\\Delta S_{N,V,E} > 0$ |\n", "| Canonical (NVT) | Helmholtz Free Energy $\\Delta A_{N,V,T} < 0$ |\n", "| Isothermal-Isobaric (NPT) | Gibbs Free Energy $\\Delta G_{N,P,T} < 0$ |\n", "\n", "\n", "In the microcanonical ensemble, meaning, energy, volume and number of particles being constant, the entropy has a simple expression that we have seen before.\n", "\n", "$$\n", "\\textrm{Entropy}: S(N,V,E) = k_B ln W\n", "$$\n", "\n", "Where W implies the number of states (the partition function if all states were equally probable). \n", "\n", "However, when instead of constant energy we keep temperature constant, that is, considering the NVT ensemble, keeping temperature, volume and number of particle constant not all states \"j\" are equally probable, so the average would be:\n", "\n", "$$\n", "\\textrm{Entropy}: S(N,V,T) = k_B \\sum_j p_j lnp_j\n", "$$\n", "\n", "and using now the canonical partition function instead of $p_j$ we get a potential that is not as elegant.\n", "\n", "$$\n", "\\textrm{Entropy}: S(N,V,T) = k_B lnQ + k_B T \\left(\\frac{\\partial ln Q}{\\partial T}\\right)_{N,V}\n", "$$\n", "\n", "with $\\Delta S_{N,V,E}>0$ as the criteria for spontaneity\n", "\n", "In the canonical ensemble, it is the free energy (Helmholtz free energy) that has a simple expression\n", "\n", "$$\n", "\\textrm{Helmholtz Free Energy}: A(N,V,T) = -k_B T ln Q\n", "$$\n", "\n", "with $\\Delta A_{V,T}<0$ as the criteria for spontaneity\n", "\n", "In the so-called isothermic isobaric ensemble where the variables held constant are N,P,T the Gibbs free energy is the thermodynamic potential with the simple expression\n", "\n", "$$\n", "\\textrm{Gibbs Free Energy}: G(N,P,T) = -k_B T ln Q\n", "$$\n", "\n", "with $\\Delta G_{P,T}<0$ as the criteria for spontaneity" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before we proceed with more mathematical expressions, it may be useful to look at real calculations of molecules and see what contributions from vibrational, translational, rotational, and electronic we should expect\n", "\n", "# Calculating thermodynamic potentials using pyscf\n", "\n", "In a separate notebook, for each of the following molecules: CO, N2, O2(use singlet!), F2, HF, HCl, CN(-), OH(-), H2.\n", "See below an example for CO.\n", "\n", "* Optimize the structure\n", "* Run a frequency analysis of the optimized structure. Save the output each in a separate file\n", " \n", "Build a function that for a given output file it reports a list of thermodynamic data such as the electronic, translational, rotational, and vibrational contributions to thermal energy, enthalpy, and free energy.\n", "\n", "Build a table with those data.\n", "\n", "Identify the most important contributions: compare what contributes the most to the internal/thermal energy, is it translational, rotational, and vibrational contributions?\n", "\n", "Compare among molecules." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "output file: m3d2_co_fq.out\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "geometric-optimize called with the following command line:\n", "/Users/pratr001/miniconda3/envs/pchem/lib/python3.12/site-packages/ipykernel_launcher.py --f=/Users/pratr001/Library/Jupyter/runtime/kernel-v3248b0861d0d945a8ae00f635235d9ec3a3ebeebc.json\n", "\n", " \u001b[91m())))))))))))))))/\u001b[0m \n", " \u001b[91m())))))))))))))))))))))))),\u001b[0m \n", " \u001b[91m*)))))))))))))))))))))))))))))))))\u001b[0m \n", " \u001b[94m#,\u001b[0m \u001b[91m()))))))))/\u001b[0m \u001b[91m.)))))))))),\u001b[0m \n", " \u001b[94m#%%%%,\u001b[0m \u001b[91m())))))\u001b[0m \u001b[91m.))))))))*\u001b[0m \n", " \u001b[94m*%%%%%%,\u001b[0m \u001b[91m))\u001b[0m \u001b[93m..\u001b[0m \u001b[91m,))))))).\u001b[0m \n", " \u001b[94m*%%%%%%,\u001b[0m \u001b[93m***************/.\u001b[0m \u001b[91m.)))))))\u001b[0m \n", " \u001b[94m#%%/\u001b[0m \u001b[94m(%%%%%%,\u001b[0m \u001b[93m/*********************.\u001b[0m \u001b[91m)))))))\u001b[0m \n", " \u001b[94m.%%%%%%#\u001b[0m \u001b[94m*%%%%%%,\u001b[0m \u001b[93m*******/,\u001b[0m \u001b[93m**********,\u001b[0m \u001b[91m.))))))\u001b[0m \n", " \u001b[94m.%%%%%%/\u001b[0m \u001b[94m*%%%%%%,\u001b[0m \u001b[93m**\u001b[0m \u001b[93m********\u001b[0m \u001b[91m.))))))\u001b[0m \n", " \u001b[94m##\u001b[0m \u001b[94m.%%%%%%/\u001b[0m \u001b[94m(%%%%%%,\u001b[0m \u001b[93m,******\u001b[0m \u001b[91m/)))))\u001b[0m \n", " \u001b[94m%%%%%%\u001b[0m \u001b[94m.%%%%%%#\u001b[0m \u001b[94m*%%%%%%,\u001b[0m \u001b[92m,/////.\u001b[0m \u001b[93m******\u001b[0m \u001b[91m))))))\u001b[0m \n", " \u001b[94m#%\u001b[0m \u001b[94m%%\u001b[0m \u001b[94m.%%%%%%/\u001b[0m \u001b[94m*%%%%%%,\u001b[0m \u001b[92m////////,\u001b[0m \u001b[93m*****/\u001b[0m \u001b[91m,)))))\u001b[0m \n", " \u001b[94m#%%\u001b[0m \u001b[94m%%%\u001b[0m \u001b[94m%%%#\u001b[0m \u001b[94m.%%%%%%/\u001b[0m \u001b[94m(%%%%%%,\u001b[0m \u001b[92m///////.\u001b[0m \u001b[93m/*****\u001b[0m \u001b[91m))))).\u001b[0m\n", " \u001b[94m#%%%%.\u001b[0m \u001b[94m%%%%%#\u001b[0m \u001b[94m/%%%%%%*\u001b[0m \u001b[94m#%%%%%%\u001b[0m \u001b[92m/////)\u001b[0m \u001b[93m******\u001b[0m \u001b[91m))))),\u001b[0m\n", " \u001b[94m#%%%%##%\u001b[0m \u001b[94m%%%#\u001b[0m \u001b[94m.%%%%%%/\u001b[0m \u001b[94m(%%%%%%,\u001b[0m \u001b[92m///////.\u001b[0m \u001b[93m/*****\u001b[0m \u001b[91m))))).\u001b[0m\n", " \u001b[94m##\u001b[0m \u001b[94m%%%\u001b[0m \u001b[94m.%%%%%%/\u001b[0m \u001b[94m*%%%%%%,\u001b[0m \u001b[92m////////.\u001b[0m \u001b[93m*****/\u001b[0m \u001b[91m,)))))\u001b[0m \n", " \u001b[94m#%%%%#\u001b[0m \u001b[94m/%%%%%%/\u001b[0m \u001b[94m(%%%%%%\u001b[0m \u001b[92m/)/)//\u001b[0m \u001b[93m******\u001b[0m \u001b[91m))))))\u001b[0m \n", " \u001b[94m##\u001b[0m \u001b[94m.%%%%%%/\u001b[0m \u001b[94m(%%%%%%,\u001b[0m \u001b[93m*******\u001b[0m \u001b[91m))))))\u001b[0m \n", " \u001b[94m.%%%%%%/\u001b[0m \u001b[94m*%%%%%%,\u001b[0m \u001b[93m**.\u001b[0m \u001b[93m/*******\u001b[0m \u001b[91m.))))))\u001b[0m \n", " \u001b[94m*%%%%%%/\u001b[0m \u001b[94m(%%%%%%\u001b[0m \u001b[93m********/*..,*/*********\u001b[0m \u001b[91m*))))))\u001b[0m \n", " \u001b[94m#%%/\u001b[0m \u001b[94m(%%%%%%,\u001b[0m \u001b[93m*********************/\u001b[0m \u001b[91m)))))))\u001b[0m \n", " \u001b[94m*%%%%%%,\u001b[0m \u001b[93m,**************/\u001b[0m \u001b[91m,))))))/\u001b[0m \n", " \u001b[94m(%%%%%%\u001b[0m \u001b[91m()\u001b[0m \u001b[91m))))))))\u001b[0m \n", " \u001b[94m#%%%%,\u001b[0m \u001b[91m())))))\u001b[0m \u001b[91m,)))))))),\u001b[0m \n", " \u001b[94m#,\u001b[0m \u001b[91m())))))))))\u001b[0m \u001b[91m,)))))))))).\u001b[0m \n", " \u001b[91m()))))))))))))))))))))))))))))))/\u001b[0m \n", " \u001b[91m())))))))))))))))))))))))).\u001b[0m \n", " \u001b[91m())))))))))))))),\u001b[0m \n", "\n", "-=# \u001b[1;94m geomeTRIC started. Version: 1.1 \u001b[0m #=-\n", "Current date and time: 2025-04-17 09:39:39\n", "#========================================================#\n", "#| \u001b[92m Arguments passed to driver run_optimizer(): \u001b[0m |#\n", "#========================================================#\n", "customengine \n", "input /var/folders/1b/pw4sp5tj3_16kjl12vvf8vt40000gq/T/tmpxbj0cqve/f3510fde-899e-48dd-b6a4-3c5e95daf0f8 \n", "logIni /Users/pratr001/miniconda3/envs/pchem/lib/python3.12/site-packages/pyscf/geomopt/log.ini \n", "----------------------------------------------------------\n", "Custom engine selected.\n", "Bonds will be generated from interatomic distances less than 1.20 times sum of covalent radii\n", "6 internal coordinates being used (instead of 6 Cartesians)\n", "Internal coordinate system (atoms numbered from 1):\n", "Distance 1-2\n", "Translation-X 1-2\n", "Translation-Y 1-2\n", "Translation-Z 1-2\n", "Rotation-A 1-2\n", "Rotation-B 1-2\n", "Rotation-C 1-2\n", " : 1\n", " : 1\n", " : 1\n", " : 1\n", " : 1\n", " : 1\n", " : 1\n", "> ===== Optimization Info: ====\n", "> Job type: Energy minimization\n", "> Maximum number of optimization cycles: 300\n", "> Initial / maximum trust radius (Angstrom): 0.100 / 0.300\n", "> Convergence Criteria:\n", "> Will converge when all 5 criteria are reached:\n", "> |Delta-E| < 1.00e-06\n", "> RMS-Grad < 3.00e-04\n", "> Max-Grad < 4.50e-04\n", "> RMS-Disp < 1.20e-03\n", "> Max-Disp < 1.80e-03\n", "> === End Optimization Info ===\n", "Step 0 : Gradient = 2.498e-01/2.498e-01 (rms/max) Energy = -112.6189499615\n", "Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 6.71931e-01\n", "Step 1 : Displace = \u001b[0m9.837e-02\u001b[0m/\u001b[0m9.837e-02\u001b[0m (rms/max) Trust = 1.000e-01 (=) Grad = \u001b[0m7.714e-02\u001b[0m/\u001b[0m7.714e-02\u001b[0m (rms/max) E (change) = -112.6652860802 (\u001b[0m-4.634e-02\u001b[0m) Quality = \u001b[0m0.998\u001b[0m\n", "Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 8.79414e-01\n", "Step 2 : Displace = \u001b[0m2.321e-02\u001b[0m/\u001b[0m2.321e-02\u001b[0m (rms/max) Trust = 1.414e-01 (\u001b[92m+\u001b[0m) Grad = \u001b[0m4.568e-02\u001b[0m/\u001b[0m4.568e-02\u001b[0m (rms/max) E (change) = -112.6663834887 (\u001b[0m-1.097e-03\u001b[0m) Quality = \u001b[0m0.324\u001b[0m\n", "Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 1.40021e+00\n", "Step 3 : Displace = \u001b[0m8.632e-03\u001b[0m/\u001b[0m8.632e-03\u001b[0m (rms/max) Trust = 1.414e-01 (=) Grad = \u001b[0m4.454e-03\u001b[0m/\u001b[0m4.454e-03\u001b[0m (rms/max) E (change) = -112.6672145698 (\u001b[0m-8.311e-04\u001b[0m) Quality = \u001b[0m1.115\u001b[0m\n", "Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 1.26368e+00\n", "Step 4 : Displace = \u001b[92m9.326e-04\u001b[0m/\u001b[92m9.326e-04\u001b[0m (rms/max) Trust = 2.000e-01 (\u001b[92m+\u001b[0m) Grad = \u001b[0m3.002e-04\u001b[0m/\u001b[92m3.002e-04\u001b[0m (rms/max) E (change) = -112.6672219101 (\u001b[0m-7.340e-06\u001b[0m) Quality = \u001b[0m0.935\u001b[0m\n", "Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 1.34887e+00\n", "Step 5 : Displace = \u001b[92m5.889e-05\u001b[0m/\u001b[92m5.889e-05\u001b[0m (rms/max) Trust = 2.828e-01 (\u001b[92m+\u001b[0m) Grad = \u001b[92m3.235e-06\u001b[0m/\u001b[92m3.235e-06\u001b[0m (rms/max) E (change) = -112.6672219430 (\u001b[92m-3.293e-08\u001b[0m) Quality = \u001b[0m0.985\u001b[0m\n", "Hessian Eigenvalues: 5.00000e-02 5.00000e-02 5.00000e-02 ... 5.00000e-02 5.00000e-02 1.34887e+00\n", "Converged! =D\n", "\n", " #==========================================================================#\n", " #| If this code has benefited your research, please support us by citing: |#\n", " #| |#\n", " #| Wang, L.-P.; Song, C.C. (2016) \"Geometry optimization made simple with |#\n", " #| translation and rotation coordinates\", J. Chem, Phys. 144, 214108. |#\n", " #| http://dx.doi.org/10.1063/1.4952956 |#\n", " #==========================================================================#\n", " Time elapsed since start of run_optimizer: 0.757 seconds\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Optimized CO structure:\n", "[[0. 0. 0.08466366]\n", " [0. 0. 1.21533634]]\n", "Rotation constant\n", "(array([ inf, 43.58846729, 43.58846729]), 'GHz')\n", "Zero-point energy\n", "(np.float64(0.002706164641888301), 'Eh')\n", "Internal energy at 0 K\n", "(np.float64(-112.6162437968851), 'Eh')\n", "Internal energy at 298.15 K\n", "(np.float64(-112.61386574425926), 'Eh')\n", "Enthalpy energy at 298.15 K\n", "(np.float64(-112.61292155970314), 'Eh')\n", "Gibbs free energy at 298.15 K\n", "(np.float64(-112.63563803079069), 'Eh')\n", "Heat capacity at 298.15 K\n", "(np.float64(8.25633742976514e-06), 'Eh/K')\n" ] } ], "source": [ "#pip install geometric\n", "\n", "from pyscf import gto\n", "from pyscf.geomopt.geometric_solver import optimize\n", "from pyscf.hessian import thermo\n", "from pyscf import lib\n", "\n", "# Convert from Bohr to Angstroms\n", "bohr_to_angstrom = lib.parameters.BOHR \n", "\n", "molCO = gto.M(\n", " atom = 'O 0 0 0; C 0 0 1.3', # in Angstrom\n", " basis = '631g',\n", " symmetry = True,\n", " charge=0, \n", " spin=0,\n", " verbose=3,\n", " output='m3d2_co_fq.out'\n", ")\n", "# Perform the Hartree-Fock calculation\n", "co_sp = molCO.RHF.run()\n", "co_opt = optimize(co_sp).run()\n", "print(\"Optimized CO structure:\")\n", "print(co_opt.atom_coords()*bohr_to_angstrom)\n", "hessian = co_sp.Hessian().kernel()\n", "\n", "# Frequency analysis\n", "freq_info = thermo.harmonic_analysis(co_sp.mol, hessian)\n", "# Thermochemistry analysis at 298.15 K and 1 atmospheric pressure\n", "thermo_info = thermo.thermo(co_sp, freq_info['freq_au'], 298.15, 101325)\n", "\n", "\n", "print('Rotation constant')\n", "print(thermo_info['rot_const'])\n", "\n", "print('Zero-point energy')\n", "print(thermo_info['ZPE' ])\n", "\n", "print('Internal energy at 0 K')\n", "print(thermo_info['E_0K' ])\n", "\n", "print('Internal energy at 298.15 K')\n", "print(thermo_info['E_tot' ])\n", "\n", "print('Enthalpy energy at 298.15 K')\n", "print(thermo_info['H_tot' ])\n", "\n", "print('Gibbs free energy at 298.15 K')\n", "print(thermo_info['G_tot' ])\n", "\n", "print('Heat capacity at 298.15 K')\n", "print(thermo_info['Cv_tot'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Questions\n", "\n", "* For the Canonical NVT ensemble use PySCF to build a table calculating the entropy, the enthalpy, and the free energy at room temperature of the molecules you already calculated in the previous M3.D3 notebook." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "pchem", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.2" } }, "nbformat": 4, "nbformat_minor": 2 }